/************************************************************************
*************************************************************************
Utility Functions
************************************************************************
************************************************************************/


/************************************************************************
Program: makebssamp
Purpose: Makes a single bootstrap dataset
************************************************************************/

cap program drop makebssamp
program makebssamp

	/******
	Step 1: Prepare the original dataset for merging
	******/
	
	//Make a household identifier that is numbered consecutively
	egen idcon = group(id)
	
	//Get the max and min value of IDs
	sum idcon
	local idmax = r(max)
	local idmin = r(min)
	
	//Get the max and min values for time
	sum t
	local tmin = r(min)
	local tmax = r(max)
	
	//Number of observations
	egen tag = tag(idcon)
	count if tag
	local nobs = r(N)
	drop tag
	
	
	//Save the source
	tempfile bssource
	save `bssource'
	
	
	/******
	Step 2: Construct framework of BS sample
	******/
	
	//Draw sample
	clear
	set obs `nobs'
	gen idcon = `idmin' + int( runiform()*(`idmax'-`idmin' + 1) )
	
	//The bootstrap household identifier
	gen id_bs = _n
	
	//Add ts
	expand `=`tmax'-`tmin'+1'
	bysort id_bs: gen t = _n + `tmin'-1
	
	tab t
	
	/******
	Step 3: Bring in data from source
	******/
	merge m:1 idcon t using `bssource'
	
	//These are households in the source that were not sampled
	drop if _m == 2
	
	//These are ts in the for which a sampled household is unobserved
	drop if _m == 1
	drop _m
	
	xtset id_bs t

end

/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end


/************************************************************************
*************************************************************************
Standard Programs
************************************************************************
************************************************************************/



/************************************************************************
Program: ar_dem
Purpose: Estimates the AR method, demeaning first
************************************************************************/
cap program drop ar_dem
program ar_dem
	preserve
	tempvar tmp
	foreach var of varlist y k l m kXk lXl mXm kXl kXm lXm {
		egen `tmp' = mean(`var')
		replace `var' = `var'-`tmp'
		drop `tmp'
	}
	
	#delimit;
	cap gmm (y - {bk}*k - {bl}*l - {bm}*m - {bkk}*kXk - {bll}*lXl - {bmm}*mXm - {bkl}*kXl - {bkm}*kXm - {blm}*lXm 
								- {rho}*(l.y - {bk}*l.k - {bl}*l.l - {bm}*l.m - {bkk}*l.kXk - {bll}*l.lXl - {bmm}*l.mXm - {bkl}*l.kXl - {bkm}*l.kXm - {blm}*l.lXm) ), 
								instruments(l.(k l m kXk-mXm) l2.m l2.k l2.kXm) 
								deriv(/bk 	= -k 	+ {rho}*l.k) 
								deriv(/bl 	= -l 	+ {rho}*l.l)
								deriv(/bm 	= -m 	+ {rho}*l.m)
								deriv(/bkk 	= -kXk 	+ {rho}*l.kXk)
								deriv(/bll 	= -lXl 	+ {rho}*l.lXl)
								deriv(/bmm 	= -mXm 	+ {rho}*l.mXm)
								deriv(/bkl 	= -kXl 	+ {rho}*l.kXl)
								deriv(/bkm 	= -kXm 	+ {rho}*l.kXm)
								deriv(/blm 	= -lXm 	+ {rho}*l.lXm)
								deriv(/rho 	= -(l.y - {bk}*l.k - {bl}*l.l - {bm}*l.m - {bkk}*l.kXk - {bll}*l.lXl - {bmm}*l.mXm - {bkl}*l.kXl - {bkm}*l.kXm - {blm}*l.lXm) ) ;
	#delimit cr
	
	restore
end




/************************************************************************
Program: ar_dem_wrap
Purpose: A wrapper for the ar_dem program (for use in bootstrapping)
************************************************************************/
cap program drop ar_dem_wrap
program ar_dem_wrap
	args spec laginst
	
	ar_dem
	if _rc == 0 {
		gen elasm_ar_dem = _b[/bm] + 2*_b[/bmm]*m + _b[/bkm]*k + _b[/blm]*l
		gen elask_ar_dem = _b[/bk] + 2*_b[/bkk]*k + _b[/bkm]*m + _b[/bkl]*l
		gen elasl_ar_dem = _b[/bl] + 2*_b[/bll]*l + _b[/blm]*m + _b[/bkl]*k
		gen ar_ar_dem	 = _b[/rho]
	}
	else {
		gen elasm_ar_dem = .
		gen elask_ar_dem = .
		gen elasl_ar_dem = .
		gen ar_ar_dem	 = .
	}



end	

/************************************************************************
Program: widstat_simple
Purpose: Estimates rho (consistently), then uses it to construct the
			Kleibergen-Paap statistic
************************************************************************/

cap program drop widstat_simple
program widstat_simple
		preserve
		sort id t
		ar_dem
		if _rc == 0 {
		local rhoest = _b[/rho]
			foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
				gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
				gen INST_`inst' = l.`inst'
			}
			
			gen yrho2 = y - `rhoest'*l.y
			
			ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id)
			local wid = e(widstat)

			restore
			gen widstat_iv = `wid'
		}
		else {
			gen widstat_iv = .
		}
	
end



/************************************************************************
Program: bspctiles_simple
Purpose: Bootstraps the percentiles of the rk-statistic
************************************************************************/

cap program drop bspctiles_simple
program bspctiles_simple
	//Number of reps, stub of the estimator (gnr or ar_dem), specification,
	//and a specifier for whether we should save or not
	args nreps savepath
	
	/******
	Step 0: Estimate true parameters
	******/
	preserve
	sort id t
	ar_dem
	local rhoest = _b[/rho]
	foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
		gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
		gen INST_`inst' = l.`inst'
	}
	
	gen yrho2 = y - `rhoest'*l.y
	
	foreach endvar of varlist ENDO2_* {
		reg `endvar' INST_* l2.m l2.k l2.kXm
		estimates store `endvar'
	}
	
	restore

	/******
	Step 1: Prepare a dataset for the bootstrapped estimates
	******/
	cap postclose bsests
	
	tempfile bsests
	postfile bsests widnull `listel' using `bsests'
	
	
	/******
	Step 2: Prepare the original dataset for bootstrapping
	******/
	tempfile original
	save `original'
	
	
	
	
	/******
	Step 3: Bootstrap
	******/
	quietly {
		forvalues b=1/`nreps' {
			preserve
				makebssamp
				drop id idcon
				sort id_bs t
				
				//Preliminary construction
				ar_dem
				
				local rhoestb = _b[/rho]
				foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
					gen ENDO2_`inst' = `inst' - `rhoestb'*l.`inst'
					gen INST_`inst' = l.`inst'
				}
				


				//Null: zero rank
				foreach endvar of varlist ENDO2_* {
					estimates restore `endvar'
					predict tmp
					replace `endvar' = `endvar' - tmp
					drop tmp
				}

				
				gen yrho2 = y - `rhoestb'*l.y
				ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id_bs)
				local widnull = e(widstat)
				

				post bsests (`widnull')
				
				/******** Display ***************/
				local disp = string(`widnull', "%9.1f")
				noisily: di "`disp' " _c
				if mod(`b',10) == 0 {
					noisily: di " "
				}
				if mod(`b',50) == 0 {
					noisily: di _col(100) "`b'"
				}
				/*******************************/
				
			restore
		}
	}
	
	postclose bsests
	
	/******
	Step 4: Compute standard errors
	******/
	
	use `bsests', clear
	
	if !mi("`savepath'") {
		save "`savepath'", replace
	}
	
	local plist 10 5 1 .5
	
	local count 0
	foreach p of numlist `plist' {
		local ++count
		egen tmp = pctile(widnull), p(`=100-`p'')
		sum tmp
		local p`count' = r(mean)
		drop tmp
	}
	
	
	use `original', clear
	
	
	foreach p of numlist 1/`count' {
		gen widstatp_`p' = `p`p''
		local pname: word `p' of `plist'
		label var widstatp_`p' "Critical value for alpha=`pname'"
	}



end



/************************************************************************
*************************************************************************
Versions that try to deal with time dummies
************************************************************************
************************************************************************/


/************************************************************************
Program: ar_dem_tdum
Purpose: Estimates the AR method, demeaning by time dummies
************************************************************************/
cap program drop ar_dem_tdum
program ar_dem_tdum
	preserve
	tempvar tmp
	foreach var of varlist y k l m kXk lXl mXm kXl kXm lXm {
		bysort t: egen `tmp' = mean(`var')
		replace `var' = `var'-`tmp'
		drop `tmp'
	}
	
	sort id t
	
	#delimit;
	cap gmm (y - {bk}*k - {bl}*l - {bm}*m - {bkk}*kXk - {bll}*lXl - {bmm}*mXm - {bkl}*kXl - {bkm}*kXm - {blm}*lXm 
								- {rho}*(l.y - {bk}*l.k - {bl}*l.l - {bm}*l.m - {bkk}*l.kXk - {bll}*l.lXl - {bmm}*l.mXm - {bkl}*l.kXl - {bkm}*l.kXm - {blm}*l.lXm) ), 
								instruments(l.(k l m kXk-mXm) l2.m l2.k l2.kXm) 
								deriv(/bk 	= -k 	+ {rho}*l.k) 
								deriv(/bl 	= -l 	+ {rho}*l.l)
								deriv(/bm 	= -m 	+ {rho}*l.m)
								deriv(/bkk 	= -kXk 	+ {rho}*l.kXk)
								deriv(/bll 	= -lXl 	+ {rho}*l.lXl)
								deriv(/bmm 	= -mXm 	+ {rho}*l.mXm)
								deriv(/bkl 	= -kXl 	+ {rho}*l.kXl)
								deriv(/bkm 	= -kXm 	+ {rho}*l.kXm)
								deriv(/blm 	= -lXm 	+ {rho}*l.lXm)
								deriv(/rho 	= -(l.y - {bk}*l.k - {bl}*l.l - {bm}*l.m - {bkk}*l.kXk - {bll}*l.lXl - {bmm}*l.mXm - {bkl}*l.kXl - {bkm}*l.kXm - {blm}*l.lXm) ) ;
	#delimit cr
	
	restore

end




/************************************************************************
Program: ar_dem_wrap
Purpose: A wrapper for the ar_dem_tdum program
************************************************************************/
cap program drop ar_dem_wrap_tdum
program ar_dem_wrap_tdum
	args spec laginst
	
	ar_dem_tdum
	if _rc == 0 {
		gen elasm_ar_dem = _b[/bm] + 2*_b[/bmm]*m + _b[/bkm]*k + _b[/blm]*l
		gen elask_ar_dem = _b[/bk] + 2*_b[/bkk]*k + _b[/bkm]*m + _b[/bkl]*l
		gen elasl_ar_dem = _b[/bl] + 2*_b[/bll]*l + _b[/blm]*m + _b[/bkl]*k
		gen ar_ar_dem	 = _b[/rho]
	}
	else {
		gen elasm_ar_dem = .
		gen elask_ar_dem = .
		gen elasl_ar_dem = .
		gen ar_ar_dem	 = .
	}



end	


/************************************************************************
Program: widstat_tdum
Purpose: Estimates rho (consistently), then uses it to construct the
			Kleibergen-Paap statistic. This version controls for time dummies
************************************************************************/

cap program drop widstat_tdum
program widstat_tdum
		preserve
		sort id t
		ar_dem_tdum
		if _rc == 0 {

			local rhoest = _b[/rho]
			foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
				gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
				gen INST_`inst' = l.`inst'
			}
			
			gen yrho2 = y - `rhoest'*l.y
			
			ivreg2 yrho2 i.t (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id)
			local wid = e(widstat)

			restore
			gen widstat_iv = `wid'
		}
		else {
			gen widstat_iv = .
		}

	
end




/************************************************************************
Program: bspctiles_tdum
Purpose: Bootstraps the percentiles of the rk-statistic
************************************************************************/

cap program drop bspctiles_tdum
program bspctiles_tdum
	//Number of reps, stub of the estimator (gnr or ar_dem), specification,
	//and a specifier for whether we should save or not
	args nreps savepath
	
	/******
	Step 0: Estimate true parameters
	******/
	preserve
	sort id t
	ar_dem_tdum
	local rhoest = _b[/rho]
	foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
		gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
		gen INST_`inst' = l.`inst'
	}
	
	gen yrho2 = y - `rhoest'*l.y
	
	foreach endvar of varlist ENDO2_* {
		reg `endvar' INST_* l2.m l2.k l2.kXm i.t
		estimates store `endvar'
	}
	
	restore

	/******
	Step 1: Prepare a dataset for the bootstrapped estimates
	******/
	cap postclose bsests
	
	tempfile bsests
	postfile bsests widnull `listel' using `bsests'
	
	
	/******
	Step 2: Prepare the original dataset for bootstrapping
	******/
	tempfile original
	save `original'
	
	
	
	
	/******
	Step 3: Bootstrap
	******/
	quietly {
		forvalues b=1/`nreps' {
			preserve
				makebssamp
				drop id idcon
				sort id_bs t
				
				//Preliminary construction
				ar_dem_tdum
				
				local rhoestb = _b[/rho]
				foreach inst of varlist k l m kXk kXl kXm lXl lXm mXm {
					gen ENDO2_`inst' = `inst' - `rhoestb'*l.`inst'
					gen INST_`inst' = l.`inst'
				}
				


				//Null: zero rank
				foreach endvar of varlist ENDO2_* {
					estimates restore `endvar'
					predict tmp
					replace `endvar' = `endvar' - tmp
					drop tmp
				}

				
				gen yrho2 = y - `rhoestb'*l.y
				ivreg2 yrho2 i.t (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id_bs)
				local widnull = e(widstat)
				

				post bsests (`widnull') 
				
				
				/******** Display ***************/
				local disp = string(`widnull', "%9.1f")
				noisily: di "`disp' " _c
				if mod(`b',10) == 0 {
					noisily: di " "
				}
				if mod(`b',50) == 0 {
					noisily: di _col(100) "`b'"
				}
				/*******************************/
				
			restore
		}
	}
	postclose bsests
	
	/******
	Step 4: Compute standard errors
	******/
	
	use `bsests', clear
	
	if !mi("`savepath'") {
		save "`savepath'", replace
	}
	
	local plist 10 5 1 .5
	
	local count 0
	foreach p of numlist `plist' {
		local ++count
		egen tmp = pctile(widnull), p(`=100-`p'')
		sum tmp
		local p`count' = r(mean)
		drop tmp
	}
	
	
	use `original', clear
	
	
	foreach p of numlist 1/`count' {
		gen widstatp_`p' = `p`p''
		local pname: word `p' of `plist'
		label var widstatp_`p' "Critical value for alpha=`pname'"
	}



end
